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Abstract — In this paper, we are interested in the classical 
problem of restoring data degraded by a convolution and 
the addition of a white Gaussian noise. The originality of 
the proposed approach is two-fold. Firstly, we formulate the 
restoration problem as a nonlinear estimation problem leading 
to the minimization of a criterion derived from Stein's unbiased 
quadratic risk estimate. Secondly, the deconvolution procedure 
is performed using any analysis and synthesis frames that can 
be overcomplete or not. New theoretical results concerning the 
calculation of the variance of the Stein's risk estimate are also 
provided in this work. Simulations carried out on natural images 
show the good performance of our method w.r.t. conventional 
wavelet-based restoration methods. 



I. Introduction 

It is well-known that, in many practical situations, one 
may consider that there are two main sources of signal/image 
degradation: a convolution often related to the bandlimited 
nature of the acquisition system and a contamination by an 
additive Gaussian noise which may be due to the electronics 
of the recording and transmission processes. For instance, the 
limited aperture of satellite cameras, the aberrations inherent to 
optical systems and mechanical vibrations create a blur effect 
in remote sensing images |2 |. A data restoration task is usually 
required to reduce these artifacts before any further processing. 
Many works have been dedicated to the deconvolution of noisy 
signals |[3l, ID, jS), |]6|. Designing suitable deconvolution 
methods is a challenging task, as inverse problems of practical 
interest are often ill-posed. Indeed, the convolution operator is 
usually non-invertible or it is ill-conditioned and its inverse is 
thus very sensitive to noise. To cope with the ill-posed nature 
of these problems, deconvolution methods often operate in a 
transform domain, the transform being expected to make the 
problem easier to solve. In pioneering works, deconvolution is 
dealt with in the frequency domain, as the Fourier transform 
provides a simple representation of filtering operations |7|. 
However, the Fourier domain has a main shortcoming: sharp 
transitions in the signal (edges for images) and other localized 
features do not have a sparse frequency representation. This 
has motivated the use of the Wavelet Transform (WT) |8l. 
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|9| and its various extensions. Thanks to the good energy 
compaction and decorrelation properties of the WT, simple 
shrinkage operations in the wavelet domain can be successfully 
applied to discard noisy coefficients flO). To take advantage 
of both transform domains, it has been suggested to combine 
frequency based deconvolution approaches with wavelet-based 
denoising methods, giving birth to a new class of restoration 
methods. The wavelet- vaguelette method proposed in |8| is 
based on an inverse filtering technique. To avoid the amplifi- 
cation of the resulting colored noise component, a shrinkage 
of the filtered wavelet coefficients is performed. The wavelet- 
vaguelette method has been refined in ifTTl by adapting the 
wavelet basis to the frequency response of the degradation 
filter However, the method is not appropriate for recovering 
signals degraded by arbitrary convolutive operators. An alter- 
native to the wavelet-vaguelette decomposition is the transform 
presented by Abramovich and Silverman |9|. Similar in the 
spirit to the wavelet-vaguelette deconvolution, a more com- 
petitive hybrid approach called Fourier- Wavelet Regularized 
Deconvolution (ForWaRD) was developed by Neelamani et 
ai: a two-stage shrinkage procedure successively operates 
in the Fourier and the WT domains, which is applicable to 
any invertible or non-invertible degradation kernel |12|. The 
optimal balance between the amount of Fourier and wavelet 
regularization is derived by optimizing an approximate version 
of the mean-squared error metric. A two-step procedure was 
also presented by Banham and Katsaggelos which employs 
a multiscale Kalman filter |13|. By following a frequency 
domain approach, band-limited Meyer's wavelets have been 
used to estimate degraded signals through an elegant wavelet 
restoration method called WaveD lfT4ll . ifTS I which is based on 
minimax arguments. In lfT6l . we have proposed an extension 
of the WaveD method to the multichannel case. 

Iterative wavelet-based thresholding methods relying on 
variational approaches for image restoration have also been 
investigated by several authors. For instance, a deconvolu- 
tion method was derived under the expectation-maximization 
framework in fH]- In flS\, the complementarity of the wavelet 
and the curvelet transforms has been exploited in a regularized 
scheme involving the total variation. In fT9\, an objective 
function including the total variation, a wavelet coefficient reg- 
ularization or a mixed regularization has been considered and a 
related projection-based algorithm was derived to compute the 
solution. More recently, the work in |20| has been extended by 
proposing a flexible convex variational framework for solving 
inverse problems in which a priori information (e.g., sparsity 
or probability distribution) is available about the representation 
of the target solution in a frame 1211 . In the same way, a 
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new class of iterative shrinkage/thresholding algorithms was 
proposed in 1*221 • Its novelty relies on the fact that the update 
equation depends on the two previous iterated values. In ||23]| . 
a fast variational deconvolution algorithm was introduced. It 
consists of minimizing a quadratic data term subject to a 
regularization on the £^-norm of the coefficients of the solution 
in a Shannon wavelet basis. Recently, in [24|, a two-step 
decoupling scheme was presented for image deblurring. It 
starts from a global linear blur compensation by a generalized 
Wiener filter Then, a nonlinear denoising is carried out by 
computing the Bayes least squares Gaussian scale mixtures 
estimate. Note also that an advanced restoration method was 
developed in ||251 . which does not operate in the wavelet 
domain. 

In the same time, much attention was paid to Stein's 
principle ||26| in order to derive estimates of the Mean- 
Square Error (MSE) in statistical problems involving an ad- 
ditive Gaussian noise. The key advantage of Stein's Unbiased 
Risk Estimate (SURE) is that it does not require a priori 
knowledge about the statistics of the unknown data, while 
yielding an expression of the MSE only depending on the 
statistics of the observed data. Hence, it avoids the difficult 
problem of the estimation of the hyperparameters of some 
prior distribution, which classically needs to be addressed in 
Bayesian approaches^ Consequently, a SURE approach can 
be applied by directly parameterizing the estimator and finding 
the optimal parameters that minimize the MSE estimate. The 
first efforts in this direction were performed in the context of 
denoising applications with the SUREShrink technique fTO|, 
||27 | and, the SUREVect estimate |28| in the case of multi- 
channel images. More recently, in addition to the estimation 
of the MSE, Luisier et al. have proposed a very appealing 
structure of the denoising function consisting of a linear 
combination of nonlinear elementary functions (the SURE- 
Linear Expansion of Threshold or SURE-LET) |29|. Notice 
that this idea was also present in some earlier works [30l. 
In this way, the optimization of the MSE estimate reduces to 
solving a set of linear equations. Several variations of the basic 
SURE-LET method were investigated: an improvement of the 
denoising performance has been achieved by accounting for 
the interscale information ||3TI and the case of color images has 
also been addressed [|32ll . Another advantage of this method is 
that it remains valid when redundant multiscale representations 
of the observations are considered, as the minimization of 
the SURE-LET estimator can be easily carried out in the 
time/space domain. A similar approach has also been adopted 
by Raphan and Simoncelli |33| for denoising in redundant 
multiresolution representations. Overcomplete representations 
have also been successfully used for multivariate shrinkage 
estimators optimized with a SURE approach operating in the 
transform domain [34 1. An alternative use of Stein's principle 
was made in ll35l for building convex constraints in image 
denoising problems. 

In |[36|, Eldar generalized Stein's principle to derive an MSE 
estimate when the noise has an exponential distribution (see 

'This does not mean that SURE approaches are superior to Bayesian 
approaches, which are quite versatile. 



also ||37]| ). In addition, she investigated the problem of the 
nonlinear estimation of deterministic parameters from a linear 
observation model in the presence of additive noise. In the 
context of deconvolution, the derived SURE was employed 
to evaluate the MSE performance of solutions to regularized 
objective functions. Another work in this research direction 
is lf38l . where the risk estimate is minimized by a Monte 
Carlo technique for denoising applications. A very recent work 
1 39 1 also proposes a recursive estimation of the risk when a 
thresholded Landweber algorithm is employed to restore data. 

In this paper, we adopt a viewpoint similar to that in f36l, 
1391 in the sense that, by using Stein's principle, we obtain an 
estimate of the MSE for a given class of estimators operating 
in deconvolution problems. The main contribution of our work 
is the derivation of the variance of the proposed quadratic risk 
estimate. These results allow us to propose a novel SURE-LET 
approach for data restoration which can exploit any discrete 
frame representation. 

The paper is organized as follows. In Section III-AI the 
required background is presented and some notations are 
introduced. The generic form of the estimator we consider for 
restoration purposes is presented in Section lTl-BI In Sectionlllll 
we provide extensions of Stein's identity which will be useful 
throughout the paper. In Section IIV-AI we show how Stein's 
principle can be employed in a restoration framework when the 
degradation system is invertible. The case of a non-invertible 
system is addressed in Section IIV-BI The expression of the 
variance of the empirical estimate of the quadratic risk is then 
derived in Section |Vl In Section [Vll two scenarii are discussed 
where the determination of the parameters minimizing the risk 
estimate takes a simplified form. The structure of the proposed 
SURE-LET deconvolution method is subsequently described 
in Section lVIll and examples of its application to wavelet-based 
image restoration are shown in Section [Villi Some concluding 
remarks are given in Section |IXl 

The notations used in the paper are summarized in Table U 

II. Problem statement 

A. Background 

We consider an unknown real-valued field whose value at 
location x e D is s(x) where D = {0, . . . , — 1} x • • • x 
{0, . . . , Dd~l} with (L>i, . . . , L>d) e {WY where N* denotes 
the set of positive integers. Here, s is a rf-dimensional digital 
random field of finite size D = Di . . . Dd with finite variance. 
Of pratical interest are the cases when d = 1 (temporal 
signals), d = 2 (images), d = 3 (volumetric data or video 
sequences) and d — A (3D+i data). 

The field is degraded by the acquisition system with (de- 
terministic) impulse response h, and it is also corrupted by 
an additive noise n, which is assumed to be independent of 
the random process s. The noise n corresponds to a random 
field which is assumed to be Gaussian with zero-mean and 
covariance field: V(x,y) S D^, E[n(x)n(y)] = 7(5x-y, where 
((5x)xez<' is the Kronecker sequence and 7 > 0. In other 
words, the noise is white. 
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TABLE I 
Notations. 



Variable 


Definition 


o 


set oi spatial (or iiequency) indices 


s 


original signal 


7, 

n 


impulse response oi the degradation tilter 


n 


additive white Gaussian noise 


r 


observed signal 


b 


rourier transiorm oi s 


TT 

ri 


Fourier transform of h 


R 


Pniiripr tran^Torm of r 

1 \JLlllCi 11 dllSHJl ill \JL 1 


F(-] 

'^y ) 


mpnn Qniinrp valiip ot thp cnal in JircriiTTipnt 

iii^ull L^UU-Ull^ \ tXlLlK^ yJL lii^ Alglldl 111 CLlgLlllll^IlL 


'^l J 


m^itlipmatipal pvnpr'tJitinTi 

liiuLllt^lliClLll^Cll t^ALfC^ldLlVJll 




Tam il V riT annlvQiQ \7PPtnrQ 
itiiiiiiy yji diiciiysij vci^ivjis 




TQTnil\? r\T c\7ntnpdc ^rpptr^vc 
id.liili y Ui sy llLllCais VCC LUi S 


) 1<£<L 


pnpffiPiPTitQ nf thp Hppnmnn'Jitinn nf q nntn ((/i/iiT^/i^r 

V^W^lll v^lK>llLi3 yjL Lll^ Li-^V^WlllL^Wol Llwll yjl. tj KJULKJ \ } 1_ ^ iL"^ Ij 




Fourier transform of (pe 




Fourier transform of ift 


e, 


estimating function applied to si 


S 


estimate oi s 


TTD 


set of frequency indices for which H is considered equal to 




set of frequency indices for which H takes significant values 


X 


threshold value in the frequency domain 


s 


projection of s onto the subspace whose 
Fourier coefficients vanish on P 




inverse Fourier transform of the projection of 

onto the subspace whose Fourier coefficients vanish on P 




index subset for the m-th subband 


A 


constant used in the Wiener-like filter 



Thus, the observation model can be expressed as follows: 

Vx £ D, r(x) = (h*s){x.)+n(x) = ^ /i(x-y),s(y)+n(x) 

yen 

(1) 

where (/i(x))xez£i is the periodic extension of (ft.(x))xGD- 
It must be pointed out that ([T]) corresponds to a periodic 
approximation of the discrete convolution (this problem can 
be alleviated by making use of zero-padding techniques Ii40il . 

ED). 

A restoration method aims at estimating s based on the 
observed data r. In this paper, a supervised approach is 
adopted by assuming that both the degradation kernel h and 
the noise variance 7 are known. 

B. Considered nonlinear estimator 

The proposed estimation procedure consists of first trans- 
forming the observed data to some other domain (through 
some analysis vectors), performing a non-linear operation on 
the so-obtained coefficients (based on an estimating function) 
with parameters that must be estimated, and finally recon- 
structing the estimated signal (through some synthesis vectors). 

More precisely, the discrete Fourier coefficients (i?(p)) 
of r are given by: 

Vp e D, i?(p) ^ r(x) exp(-27rzx^D"^p) (2) 

where D — Diag(I?i, . . . , Dd). In the frequency domain, (HJ 
becomes: 

R{p)^U{p) + N{p), where U (p) = H (p) S (p) (3) 



and, the coefficients S{p) and A^(p) are obtained by expres- 
sions similar to (|2]i. 

Let {(pi)i<i<L be a family of L € N* analysis vectors 
of M'°i^ -^^''7Thus, every signal r of R^i^ '^-^'' can be 
decomposed as: 

V^e {!,...,£}, -(r,^,) =^r(x)^,(x), (4) 

xeD 

the operator (•, •) designating the Euclidean inner product of 
]gDix x_Dd According to Plancherel's formula, the coeffi- 
cients of the decomposition of r onto this family are given 
by 

V^G r, = -l^i?(p)($,(p))*, (5) 

where ^e{p) is a discrete Fourier coefficient of (pg and (•)* 
denotes the complex conjugation. Let us now define, for every 
£ e {!...., L}, an estimating function 9^ : M ^ K (the choice 
of this function wiU be discussed in Section IVII-BI ). so that 

S£ = e,(r,). (6) 

We will use as an estimate of s(x), 

L 

s(x) = ^S£^£(x) (7) 

(=1 

where {(fii)i<e<L is a family of synthesis vectors of 
jjDix - xDd Equivalently, the estimate of S is given by: 

L 

^(p) = ^?,$,(p) (8) 

£=1 
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where $£(p) is a discrete Fourier coefficient of (^i. It must 
be pointed out that our formulation is quite general. Different 
analysis/synthesis families can be used. These families may 
be overcomplete (which implies that L > D) or not. 

III. Stein-like identities 

Stein's principle will play a central role in the evaluation 
of the mean square estimation error of the proposed estimator. 
We first recall the standard form of Stein's principle: 

Proposition 1. |26| Let Q:M.~^Rbea continuous, almost 
everywhere dijferentiable function. Let rj be a real-valued zero- 
mean Gaussian random variable with variance and v be 
a real-valued random variable which is independent of rj. Let 
p — V + 7] and assume that 



SI. 



0, 



Vr e M, lim|(;|^oo B(t + () exp I 

E[(e(p))^] < oo and E[|e'(p)|] < oo where Q' is the 
derivative of Q. 



Then, 



(9) 



We now derive extended forms of the above formula (see 
Appendix |A]i which will be useful in the remainder of this 
paper: 

Proposition 2. Let 9j : R ^ R with i G {1,2} be continuous, 
almost everywhere differentiable functions. Let {rji, »72, '72) 
be a real-valued zero-mean Gaussian vector and (^1,^2) 
be a real-valued random vector which is independent of 
iVii V2i'nij ^2)- Pi = Vi -\- rji where i € {1,2} and assume 
that 

A)- 

2a^ I 



Vr e M, limi 



CI- 



e,(r+C)C'exp(- 



(i) Vae 
0, 

(ii) E[|e,(pOI'] <oo, 

(iii) E[|0J(/9i)p] < 00 where Q[ is the derivative of Qi. 
Then, 



E[ei(pi)^yi] 

E[9l(/C»l)^l7y2] 

E.[Qi{pi)Q2{p2)mm] 



(10) 



(11) 
2E[e'i(pi)] 

(12) 



:E[eUpi)]E[77i^yi] 

X E[77i^2]E[?72?7i]- 

+ E[ei(pi)e^(p2)^?i]E[r;2^2] 

+ E[eUpi)e^(p2)](E[ryi?y2]E[772m] 

- E[77i^i]E[772?72])- (13) 

Note that Proposition |2] obviously is applicable when 
{vi,V2) is deterministic. 

IV. Use OF Stein's PRINCIPLE 

A. Case of invertible degradation systems 

In this section, we come back to the deconvolution problem 
and develop an unbiased estimate of the quadratic risk: 



(14) 



which will be useful to optimize a parametric form of the esti- 
mator from the observed data. For this purpose, the following 
assumption is made: 

Assumption 1. 

(i) The degradation filter is such that, for every p G D, 
i/(p)^0. 

(ii) For every £ in {1, . . . , L}, Qe is a continuous, almost 
everywhere differentiable function such that 

(a) Va e M*,Vt e R, 

lim e,(r + C)C'exp(-f^) = 0, 

(b) E[\Q(,{r(,)\^] < 00 and E[|e^(r<,)|3] < 00 where 
0^ is the derivative of Qg. 

Under this assumption, the degradation model can be re- 
expressed as s(x) = r(x) — n(x) where r and n are the fields 
whose discrete Fourier coefficients are 



(15) 



i7(p)' Hip)- 

Thus, since the noise has been assumed spatially white, it is 
easy to show that 



V(p,p') e 



E[iV(p)(7V(p'))' 



-yD 



Sp-p' (16) 



Hip) 

E[iV(p)(Ar(p'))*] -^^Vp' 



(17) 

and E[7V(p)(S'(p'))*] — 0. The latter relation shows that n 
and s are uncorrelated fields. 

We are now able to state the following result (see Appendix 
D: 

Proposition 3. The mean square error on each frequency 
component is such that, for every p g D, 

E[|^(p)-5(p)p] = E[|^(p)-i?(p)p]- 



I^(P) 



+ 2,X:E[eK..)]Re{^M|M:} (18) 
i=l ^P-^ 

and, the global mean square estimation error can be expressed 
as 

E[£:(s-s)] = E[£(s-?)] + A (19) 

L 

where (7£)i<^<l is the real-valued cross-correlation sequence 
defined by: for all £ G {I, ... , L}, 



1 V- *£(P)(*^(P))' 



D ^ H(p) 

pea 



(21) 



B. Case of non-invertible degradation systems 

Assumption[]Ui)|expresses the fact that the degradation filter 
is invertible. Let us now examine how this assumption can be 
relaxed. 
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We denote by P the set of indices for which the frequency 
response H vanishes: 



= {p e D I i?(p) - 0}. 



(22) 



It is then clear that the components of 5(p) with p e P, 
are unobservable. The observable part of the signal s thus 
corresponds to the projection s = n(s) of s onto the subspace 
of X ■ ■ X Dd j.jjg fields whose discrete Fourier coefficients 
vanish on P. In the Fourier domain, the projector 11 is therefore 
defined by 



Vp e 



^(p) ifp^: 
if p e 







(23) 



In this context, it is judicious to restrict the summation in (|5]l 
to Q = ]D)\P so as to limit the influence of the noise present in 
the unobservable part of s. This leads to the following modified 
expression of the coefficients r^: 



-l^i?(p)(<i>,(p))^ 

pGQ 



(24) 



where r = n(r). The second step in the estimation procedure 
(Eq. ^) is kept unchanged. For the last step, we impose the 
following structure to the estimator: 



(25) 



where Lp^ — n((^^). We will also replace Assumption [lUi)| by 
the following less restrictive one: 

Assumption 2. The set Q is nonempty. 



Under this condition and Assumption \ ^\\)\ an extended 
form of Proposition [3] is the following: 

Proposition 4. The mean square error on each frequency 
component is given, for every p € Q, /^y (II 81 ). The global 
mean square estimation error can be expressed as 



£[£(s- s)] = £[£{s - s)] + £[£{s-r)] + A (26) 

L 

^ = ^(2E - E l^(P)r')- (27) 

1=1 p6Q 

Hereabove, r_ denotes the 2D field with Fourier coefficients 



R{p) = < ^(p) 

[O otherwise, 



(28) 



and, the real-valued cross-correlation sequence {'jg)i<i<L 
becomes: _ 

_ 1 *^(P)(*^(P))* 



^(P) 



(29) 



Proof: The proof that ( fTSl ) holds for every p G Q is 
identical to that in Proposition [3] The global MSE can be 
decomposed as the sum of the errors on its unobservable 
and observable parts, respectively. Using the orthogonality 
property for the projection operator 11, the corresponding 
quadratic risk is given by: £{s — s) £{s — s) + £{s ~ s). 



It remains now to express the mean square estimation error 
E[£(J— s)] on the observable part. This is done quite similarly 
to the end of the proof of Proposition [3] ■ 

Remark 1. 

(i) Assume that the functions s and h share the same 
frequency band in the sense that, for all p ^ Q, 
5(p) = 0. Let us also assume that, for every p G Q, 
H{p) ~ 1. This typically corresponds to a denoising 
problems for a signal with frequency band Q. Then, 
since s = s and r = r. i26i becomes 



+ -^(2)jE[eKr,)](^,,^;-car 



L 



(30) 



where card(Q) denotes the cardinality of Q. In the 
case when d = 2 (images) and Q = D, the resulting 
expression is identical to the one which has been derived 
in [29] for denoising problems. 
(ii) Proposition |4] remains valid for more general choices 
of the set P than ( 122b . In particular, (126b and ( 127b are 
unchanged if 



-{pe 



l^^(p)l < x} 



(31) 



where x ^ 0, provided that the complementary set Q 
satisfies Assumption \2\ 
(iii) It is possible to give an alternative proof of ( I26I )-( I27I ) 
by applying Proposition 1 in A36i . 

V. Empirical estimation of the risk 

Under the assumptions of the previous section, we are now 
interested in the estimation of the "observable" part of the risk 
in (|26] |. that is £o = £(3 — s), from the observed field r. As 
shown by Proposition |4l an unbiased estimator of £0 is 



£o=£{s-r) + A 



(32) 



where 



^ = ^(2Eq^(^^)7£ - E i^(p)r')- (33) 

i=i peQ 

We will study in more detail the statistical behaviour of this 
estimator by considering the difference: 

- - ^ E (^(^) - ^(^)) ^(^) - ^® - (34) 

More precisely, by making use of Proposition |2] the variance 
of this term can be derived (see Appendix 0. 

Proposition 5. The variance of the estimate of the observable 
part of the quadratic risk is given by 

yar[£o~£o] = ^E[£{sh -th)] 



E E minmn)]%,^.,-% E mm 



pet; 



(35) 
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where th is the field with discrete Fourier coefficients given 

by 

^.(P)4|iy '^^^^ (36) 
I otherwise, 

s'h is similarly defined from 's and, 

y{£,i) e{l,...,L} , 7^,, ==7^2^ 



pet 



^(P) 



(37) 



Remark 2. 

(i) Eq. ( I35l l suggests that caution should be taken in re- 
lying on the unbiased risk estimate when \H{p)\ takes 
small values. Indeed, the terms in the expression of the 
variance involve divisions by H(jp) and may therefore 
become of high magnitude, in this case. 

(ii) An alternative statement of Proposition |5] is to say that 



47, 



L L 



£|2 



pet; 



is an unbiased estimate of Var[£o — £o]- 

VI. Case study 

It is important to emphasize that the proposed restoration 
framework presents various degrees of freedom. Firstly, it 
is possible to choose redundant or non redundant analy- 
sis/synthesis families. In Section IVI-AI we will show that 
in the case of orthonormal synthesis families, the estimator 
design can be split into several simpler optimization proce- 
dures. Secondly, any structure of the estimator can be virtu- 
ally considered. Of particular interest are restoration methods 
involving Linear Expansion of Threshold (LET) functions, 
which are investigated in Section IVI-B I As akeady mentioned, 
the latter estimators have been successfully used in denoising 
problems [i29i . 

A. Use of orthonormal synthesis families 

We now examine the case when {'fj)i<g<L is an orthonor- 
mal basis of n(R^i^ •■■^^<') (thus, L = card(Q)). This arises, 
in particular, when {(pe)i<e<L is an orthonormal basis of 
jjDi x - xDd ^jjj j.jjg degradation system is invertible (Q — D). 
Then, due to the orthogonality of the functions {(f^)i<i<L, 
the unbiased estimate of the risk in ( |26] | can be rewritten as 
£{s — s)+£o = S{s — s) + D^^ J2e=ii^^ ~Le)'^ + where 
Zi = {l, 'fi)- Thanks to (|33] |, the observable part of the risk 
estimate can be expressed as 



£) - -'' D ^ " " D 

1=1 1=1 pec 



where (7^)i<^<i is given by (|29] l. 



(38) 



Let us now assume that the coefficients {ri)i<i<L are 
classified according to M e N* distinct nonempty index 
subsets IK,„, TOG {!,..., A/}. We have then L = X]m=i 
where, for every m £ {!,.... M}, Km — card(]K,„). For 
instance, for a wavelet decomposition, these subsets may 
correspond to the subbands associated with the different res- 
olution levels, orientations,... In addition, consider that, for 
every m e {!,..., A/}, the estimating functions {Qe)^^^ 
belong to a given class of parametric functions and they are 
characterized by a vector parameter a„j. The same estimating 
function is thus employed for a given subset K„i of indices. 
Then, it can be noticed that the criterion to be minimized in 
(|38] | is the sum of M partial MSEs corresponding to each 
subset IK„j. Consequently, we can separately adjust the vector 
am, for every me {!,..., A/}, so as to minimize 



E (©H^") 



27 E 0^(^^)7£- (39) 



B. Example of LET functions 

As in the previous section, we assume that the coefficients 
(?'^)i<£<L as defined in (l24l i are classified according to 
M eW distinct index subsets Km, m e {1, . . . , M}. Within 
each class K^, a LET estimating function is built from a 
linear combination of Im G N* given functions fm.i- M ^ M 
applied to r^. So, for every me {!,..., Af } and £ e K^, the 
estimator takes the form: 



(40) 



where {am,i)i<i<i^^ are scalar real-valued weighting factors. 
We deduce from dZST l that the estimate can be expressed as 



M Irr 



where 



^m.zW " E fmA'rt)'£f{^)- 



(41) 



(42) 



Then, the problem of optimizing the estimator boils down to 
the determination of the weights am,i which minimize the 
unbiased risk estimate. According to (|32] | and (l33T l. this is 
equivalent to minimize f (s-£) + ^ I]m=i S^eK^ ^'lij^dlf.^ 
where (7^)]^<^<^ is given by ( [29l ). From (ITTT l. it can be 
deduced thatThrs amounts to minimizing: 

E E ""0''o E E°™''^<^™o,^o'^m,^^ 
mo — 1 20 — 1 m— 1 i— 1 

-2 E«™o..o(^^„,,„,i) 

mo = l io = l 

+ E E E "'"«.»0'^™0,2o(^^)7^- 

mo — 1 i£¥^jnf. 20 — 1 



7 



This minimization can be easily shown to yield the following 
set of linear equations: 

Vmo e {l,...,Af},Vio G {!,...,/„„}, 

' Clm..i 



Z — Z — j^'—mo,to —n 
ni — 1 i—1 

(^™o.o'£)-^ E /™o,.o(-^)7.. (43) 



VII. Parameter choice 



A. Choice of analysis/synthesis functions 

Using the same notations as in Section IVll let {K„i,l < 
m < M} be a partition of {1,...,L}. Consider now a 
frame of M.d^x--xd^^ ((V'm.kJ^eK^) i<,„<^^, where, for ev- 
ery m g {1, . . . ,M}, ■(/'m^o is some field in R-Dix - xDd ^j^^^ 
for every £ e K™, ipm^i denotes its -periodically shifted 
version where is some shift value in D. Notice that, by 
appropriately choosing the sets (K„i)i<„i<m, any frame of 
x - xDd can be written under this form but that it is mostly 
useful to describe periodic wavelet bases, wavelet packets ll43l . 
mirror wavelet bases [HI, redundant/undecimated wavelet 
representations as well as related frames [44 1, ||45l , ||451 , ||47l . 
For example, for a classical ID periodic wavelet basis, M —1 
represents the number of resolution levels and, for every £ in 
subband K.,„ at resolution level me {1, . . . , M — 1}, the shift 
parameter k£ is a multiple of 2™ (Km being here the index 
subset related to the approximation subband). 

A possible choice for the analysis family {(pe)i<£<L is then 
obtained by setting 

Vm e {1, . . . , M}, y£ e K,„, Vp e B, 

$£(p) = G(p)^'™,k,(p) 

= exp(-2^ik/D-ip)G(p)*„,.o(p) 

(44) 

where G{p) typically corresponds to the frequency response 
of an "inverse" of the degradation filter. It can be noticed 
that a similar choice is made in the WaveD estimator 1 14] by 
setting, for every p e Q, G(p) ~ 1/ (if(p)) (starting from a 
dyadic Meyer wavelet basis). By analogy with Wiener filtering 
techniques, a more general form for the frequency response of 
this filter can be chosen: 



G(p) = 



H{P) 



(45) 



|i/(p)P + A 

where A > 0. Note that, due to (l44l l. the computation of 
the coefficients (r£)i<£<L amounts to the computation of the 
frame coefficients: 

Vm e {1, . . . , M},W e K„, n = {r, V-m^k.) (46) 
where f is the field with discrete Fourier coefficients 



^(p)^ f(G(p))*i?(p) ifpeQ 
I otherwise. 



(47) 



Concerning the associated synthesis family 
{^i)i<e<L, we simply choose the dual synthesis 



frame of ((V'm.kf )fGK„) i<„<j,^which, with a slight 
abuse of notation, will l)e assumed of the form: 
{{'?m,ke)i&K^) i<„i<M where, for every £ e K,,,., ip„i,ke 
denotes the kf-periodically shifted version of ^m o- So, 
basically the restoration method can be summarized by Fig. 

m 



With these choices, it can be deduced from 



that 



Vme {i,...,M},y£e 



-E 



peQ 



1'm.o(p)(^m,o(p))'' 

I^(P)P + A 



(48) 



This shows that only M values of 7^ need to be computed 
(instead of L). Similarly, simplified forms of the constants 
(7^ j)i<£^i<L and (k^ )i<£<L as defined by dJTl ) and ( 1124b can 
be easily obtained. 

B. Choice of estimating functions 

We will employ LET estimating functions due to the sim- 
plicity of their optimization, as explained in Section IVI-BI 
More precisely, the following two possible forms will be 
investigated in this work: 

• nonlinear estimating function in ||29l : we set /,„ = 2, 
take for /„i_i the identity function and choose 



P (49) 



VpeR,/™.2(p)= l-cxp - 



where lo e]0,cx)[ and (t„j is the standard deviation 
of («£)£gK,„- According to (|120| i and (|44] |. we have, 
for any /"e K„, = 7^"' EpeQ l'i'£(p)l' = 
7i)-^EpeQ|G(p)Pl^™,o(p)p. 
• nonlinear estimating function in ll30l : again, we set /,„ = 
2, and take for f„i.i the identity function but, we choose: 

Vp e M, 

I,na{p) = ( tanh (^^) - tanh (P^)) p 

(50) 

where G]0;Oop and dm is defined as for the 

previous estimating function. 

VIII. Experimental results 
A. Simulation context 

In our experiments, the test data set contains six 8-bit images 
of size 512 x 512 which are displayed in Fig. |2] Different 
convolutions have been applied: (i) 5 x 5 and 7x7 uniform 
blurs, (ii) Gaussian blur with standard deviation equal to 
2, (iii) cosine blur defined by: V(pi,p2) G {0, . . . , I?i — 1} x 
{0, . . . ,i?2 - 1}, li{vuV2) = Hi{pi)H2{p2) where 



Vze{l,2}, 



cos 



7r(pi - FcDi) 



V (1 



(1 - 2F,)D, 
(i/,(A-K))' 



if < < F,D, 
if FcD, <p,< A/2 
otherwise 

(51) 
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Fig. 1. Restoration metliod. 



with Fc £ [0, 1/2), (iv) Dirac (the restoration problem then re- 
duces to a denoising problem) and, realizations of a zero-mean 
white Gaussian noise have been added to the blurred images. 
The noise variance 7 is chosen so that the averaged blurred 
signal to noise ratio BSNR reaches a given target value, where 
BSNR= lOlogio [\\h*s ||2 /{D'y)y The performance of a 
restoration method is measured by the averaged Signal to 
Noise Ratio: SNR= lOlogm (e[s2]/E[(s - s)^]^ where E 
denotes the spatial average operator. In our simulations, we 
have chosen the set P as given by dsTl i where the threshold 
value X is automatically adjusted so as to secure a reliable 
estimation of the risk while maximizing the size of the set Q. 
In practice, x h^s been set, through a dichotomic search, to 
the smallest positive value such that £0 > lOVKnax, where 
Vmax is an upper bound of Var[£o — £o]- This bound has 
been derived from ( |35] ) under some simplifying assumptions 
aiming at facilitating its computation. In an empirical manner, 
the parameter A in ( |45] ) has been chosen proportional to the 
ratio of the noise variance to the variance of the blurred 
image, by taking A = 37/(E[r^] — (E[r])^ — 7). The other 
parameters of the method have been set to cj = 3 in ( |49l ) and 
= (3.5,2.25) in §0^. 
To validate our approach, we have made comparisons with 
state-of-the-art wavelet-based restoration methods and some 
other restoration approaches. For all these methods, symlet- 
8 wavelet decompositions performed over 4 resolution levels 
have been used |42|. The first approach is the ForWaRD 
methocjl which employs a translation invariant wavelet repre- 
sentation 1481 . i49l . The ForWaRD estimator has been applied 
with an optimized value of the regularization parameter The 
same translation invariant wavelet decomposition is used for 
the proposed SURE-based method. The second method we 
have tested is the TwISlH algorithm |22| considering a total 
variation penalization term. The third approach is the varia- 
tional method in li2Tl Section 6] (which extends the method 

Matlab toolbox can be downloaded from 

http://www.dsp.rice.edu/software/ward.shtml. 

Matlab toolbox can be downloaded from 
http://www.lx.it.pt/~bioucas/code.htm. 



in f20l) where we use a tight wavelet frame consisting of 
the union of four shifted orthonormal wavelet decompositions. 
The shift parameters are (0,0), (1,0), (0,1) and (1,1). We 
have also included in our comparisons the results obtained with 
the classical Wiener filter and with a least squares optimization 
approach using a Laplacian regularization operator^ 

B. Numerical results 

Table provides the values of the SNR achieved by the 
different considered techniques for several values of the BSNR 
and a given form of blur (uniform 5 x 5) on the six test 
images. All the provided quantitative results are median values 
computed over 10 noise realizations. It can be observed that, 
whatever the considered image is, SURE-based restoration 
methods generally lead to significant gains w.rt. the other 
approaches, especially for low BSNRs. Furthermore, the two 
kinds of nonlinear estimating function which have been eval- 
uated lead to almost identical results. It can also be noticed 
that the ForWaRD and TwIST methods perform quite well 
in terms of MSE for high BSNR. However, by examining 
more carefully the restored images, it can be seen that these 
methods may better recover uniform areas, at the expense of 
a loss of some detail information which is better preserved by 
the considered SURE-based method. This behaviour is visible 
on Fig. [3] where the proposed approach allows us to better 
recover Barbara's stripe trouser. 

Table |III] provides the SNRs obtained with the different 
techniques for several values of the BSNR and various blurs 
on Tunis image (see Fig.|2](e)). The reported results allow us to 
confirm the good performance of SURE-based methods. The 
lower performance of the wavelet-based variational approach 
may be related to the fact that it requires the estimation of 
the hyperparameters of the prior distribution of the wavelet 
coefficients. This estimation has been performed by a max- 
imum likelihood approach which is suboptimal in terms of 
mean square restoration error. The results at the bottom- 
right of Table |III] are in agreement with those in ||29l , ||33l 

''We use the implementations of these methods provided in the Matlab 
Image Processing Toolbox, assuming that the noise level is known. 
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showing the outperformance of LET estimators for denoising 
problems. The poorer results obtained with ForWaRD in this 
case indicate that this method is tailored for deconvolution 
problems. 

In the previous experiments, for all the considered methods, 
the noise variance 7 was assumed to be known. Table |IV] 
gives the SNR values obtained with the different techniques 
for several noise levels and various blurs on Tunis image, 
when the noise variance is estimated via the classical median 
absolute deviation (MAD) wavelet estimator ll50l p. 447]. One 
can observe that the results are close to the case when the 
noise variance is known, except when the problem reduces 
to a denoising problem associated with a high BSNR. In this 
case indeed, the MAD estimator does not provide a precise 
estimation of the noise variance. However, the restoration 
results are still satisfactory. 

IX. Conclusions 

In this paper, we have addressed the problem of recovering 
data degraded by a convolution and the addition of a white 
Gaussian noise. We have adopted a hybrid approach that com- 
bines frequency and multiscale analyses. By formulating the 
underlying deconvolution problem as a nonlinear estimation 
problem, we have shown that the involved criterion to be 
optimized can be deduced from Stein's unbiased quadratic risk 
estimate. In this context, attention must be paid to the variance 
of the risk estimate. The expression of this variance has been 
derived in this paper 

The flexibility of the proposed recovery approach must 
be emphasized. Redundant or non-redundant data represen- 
tations can be employed as well as various combinations of 
linear/nonlinear estimates. Based on a specific choice of the 
wavelet representation and particular forms of the estimator 
structure, experiments have been conducted on a set of images, 
illustrating the good performance of the proposed approach. 
In our future work, we plan to further improve this restoration 
method by considering more sophisticated forms of the estima- 
tor, for example, by taking into account multiscale or spatial 
dependencies as proposed in [32J, |34J for denoising problems. 
Furthermore, it seems interesting to extend this work to the 
case of multicomponent data by accounting for the cross- 
channel correlations. 

Appendix A 
Proof of Proposition[2] 

We first notice that Assumption |(ii)| is a sufficient condition 
for the existence of the left hand-side terms of (fT0t-(fT3]l since. 



by Holder's inequality, 

mi{Pl)m\] <E[|ei(pi)|3]l/3E[|?/i|3/2]2/3 (52) 
E[|ei(pi)^yi^72|]<E[|ei(pi)|3]l/3E[|?^l^,|3/2]2/3 

(53) 

E[|ei(pi)?/i|^]<E[|ei(pi)|3]l/3E[|?/i|3/2|5y2|3]2/3 

(54) 

E[|ei(pi)ei(p2)m^2|]<E[|ei(pi)|3]i/3E[|ei(p2)|^]^/^ 

X E[\m?j2\'']'/'. (55) 
We can decompose rji with i E {1,2} as follows: 

Vi = 0^771 + f]i (56) 
where a.i is the mean-square prediction coefficient given by 

a^ai = E[riirji] (57) 

with (7^ = E[?7^] and, fji is the associated zero-mean prediction 
error which is independent of vi and 771 lf| We deduce that 

E[ei(pi)^i] =aiE[ei(pi)r;i]. (58) 

We can invoke Stein's principle to express E[0i(pi)77i], pro- 
vided that the assumptions in Proposition [T] are satisfied. To 
check these assumptions, we remark that, for every t G M, 
when 1^1 is large enough, |9i(r + C)| exp (— ^) < |9i(r + 
C)|C^exp ( — ^3"), which, owing to Assumption [(1)] implies 
that lim|^|_^oo + C) exp ( — ^3-) = 0. In addition, 

from Jensen's inequality and Assumption |(ii)[ E[|8']^(pi)|] < 
E[|6i(pi)p]^/^ < 00. Consequently, Q combined with ( l57b 
can be applied to simplify ( fSSl ). so allowing us to obtain ( fTOb . 
Let us next prove ( fTTT ). From ( l56b , we get: 

E[6i(pi)^i7y2] =aia2E[ei(pi)?7i] + aiE[ei(/9i)77i]E[772] 

+ a2E[ei(pi)77i]E[r7i] + E[ei(pi)]E[77i772] 
=aia2E[ei(pi)772] + E[ei(pi)]E[77i772] (59) 

where we have used in the first equality the fact that (771 ,772) is 
independent of (771 , vi ) and, in the second one, that it is zero- 
mean. Then, by making use of the orthogonality relation: 

E[^i^2] = aia2(T^ + E[7)i772] (60) 

we have 

E[ei(pi)^i^2] =aia2(E[ei(pi)772] - c7^E[ei{pi)]) 

+ E[ei(pi)]E[^i^2]. (61) 

'Recall that , »7i , f?2 ) is zero-mean Gaussian. 
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TABLE II 

Restoration results for a 5 x 5 uniform blur: initial SNR (SNRJ) and SNR obtained with our approach using the nonlinear 
function in ^49) (snr_b), our approach using the nonlinear function in jso) (snr_s), forward (snr_f), twist (snr_t), the 
wavelet-based variational approach (snr_v), the wiener filter (snr_w) and the regularized quadratic method (snr_r). 
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20.24 




SNR_b 


17.02 
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19.75 




SNR_b 


18.63 


19.73 


20.75 


21.72 


22.66 




SNR_s 


16.99 


17.52 


18.06 


18.77 


19.63 




SNR_s 


18.62 


19.73 


20.73 


21.70 


22.66 


Barbara 


SNRJ 


16.14 


17.04 


17.62 


18.56 


19.49 
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SNRJ 


16.57 
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19.99 
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18.10 




SNR_w 


15.80 


19.07 


20.37 


20.79 
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SNR_w 


12.60 


14.70 


15.53 


15.81 


15.90 




SNR_w 


18.01 


23.13 


25.48 


26.28 


26.53 




SNR_r 


13.25 


14.43 


15.45 


15.78 


15.90 




SNR_r 


23.65 


24.60 


25.42 


26.12 


26.44 




In addition, by integration by parts, the conditional expectation can be reexpressed as 
w.r.t. vi given by r\n, / \ '2 \ i 



Vt, EiOMrjl \ VI ^ t] 



= ( lim 8i(T + C)Cexpi 



lim ei(r + C)Cexp(-^) 



1 c 

—j_j,ir^oe^M-y< (62) + r(e,,. + o + e;(. + cK)exp(-^)<ic). 



(63) 
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TABLE III 

Tunis image restoration for various blurs. 
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ffh = 2 
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7x7 


SNR_f 
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SNR_t 
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SNR_t 
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SNR_b 


19.97 


22.30 


25.05 


28.25 


31.97 
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18.52 
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SNR_t 


19.82 
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31.00 
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SNR_v 


18.10 


21.02 
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30.18 




SNR_w 


12.63 


17.50 


21.98 


25.61 


27.89 




SNR_w 


10.42 


15.13 


20.04 


25.01 


30.00 




SNR_r 


19.23 


21.03 


23.00 


25.03 


27.04 




SNR_r 


19.39 


21.37 


23.76 


26.58 


29.91 



TABLE IV 

Tunis image restoration for various blurs when the noise variance is estimated with the MAD estimator. 



Blur 


BSNR 


10 


15 


20 


25 


30 


Blur 


BSNR 


10 


15 


20 


25 


30 




SNR_i 


9.662 


13.86 


17.01 


18.80 


19.56 




SNR_i 


9.577 


13.64 


16.56 


18.13 


18.77 


Gaussian 


SNR_b 


18.16 


19.12 


20.00 


20.82 


21.57 


Unifomi 


SNR_b 


18.05 


18.87 


19.57 


20.30 


21.15 


an = 2 


SNR_s 


18.15 


19.11 


19.99 


20.81 


21.56 


7x7 


SNR_s 


18.05 


18.86 


19.57 


20.30 


21.15 




SNR_i 


9.971 


14.87 


19.57 


23.74 


26.83 




SNR_i 


10.00 


15.00 


20.00 


25.00 


30.00 


Cosine 


SNR_b 


19.82 


21.87 


24.05 


26.15 


27.99 


Dirac 


SNR_b 


19.97 


22.28 


24.96 


27.87 


30.79 


Fc = 3/32 


SNR_s 


19.81 


21.85 


24.02 


26.12 


27.97 




SNR_s 


19.96 


22.26 


24.92 


27.83 


30.77 



The existence of the latter integral is secured for almost every 
value T that can be taken by vi, thanks to Assumptions |(ii)| 
and |(iii)| and, the fact that, if /x denotes the probability measure 

of Vi, 




=E[|ei(pi) + e;(pi)ryi|] 

<E[|ei(pi)|] + E[|el(pi)r?i|] 

<E[\eM\Y^' + E[|el(pi)P]^/'E[hip/2]2/3 < oo. 

(64) 

Since, for every r e M, when |C| is large enough, |8i(r + 
OCI exp {-£2) < |ei(T + C)|C' exp ( - ^), Assumption 
^ implies that lim|^|^oo 0i(t + OCexp ( - ^) = 0. By 
using this property, we deduce from ( l63T l that £[81(^1)77^ | 
vi] = o-2(E[ei(wi+77i) I vi] + E[e[{vi+r]i)vi I ^i]), which 
yields 

E[ei(pi)r;2] = a2(E[ei(pi)] + E[e'Mm])- (65) 

By inserting this equation in $6T[ . we find that 
E[ei{pi)?ji?j2] = aia2a^E[Q[{pi)rji] + E[ei(pi)]E[^i^y2]. 
Formula ( fTTT i straightforwardly follows by noticing that, 
according to (|56] l and (|57] i, E[6'j(/9i)772]E[77i?yi] = 
aia2(J^E[Q'i{pi)T]i]. Consider now (fT2T i. By using (|56] l 
and the independence between (771,7/2) and (r/i,ui), we can 



write 

Ei&Mmvi] 

^aialE[ei{pi)r]l] + 2aia2E[Qi{pi)Tjl]E[fi2] 
+ aiE[Qi{pi)m]E[fii] + alE[eMTjf]E[fii] 
+ 2a2E[ei(pi)7?i]E[77i772] + E[ei(pi)]E[77i?7^] 

=aialE[ei{pi)vf] + E[ei(pi)77i](aiE[772] + 2a2E[mfi2]) 

(66) 

where the latter equality stems from the symmetry of the 
probability distribution of (771,772)- Taking into account the 
relation E[rj^] ~ a^cr^ + £[77!] and ( I6OI 1, we get 

E[ei(pi)^7i^7|] =aia2 (£[61(^1)77?] - 3a'E[eMm]) 
+ E[ei(pi)77i](aiE[^] + 2a2E[^?i^y2]). 

(67) 

Let us now focus our attention on E[Oi{pi)r]l]. Since As- 
sumptions |(ii)| and |(iii)| imply that 

E[|2ei(pi)77i+e;(pi)7y2|] < 2E[\eM\']'/'E[\m\'^'r^' 

+ E[|e'i(pi)|3]l/3E[|^l|3]2/3<^ (68) 

and Assumption [(1)] holds, we can proceed by integration by 
parts, similarly to the proof of ( |65] l to show that 

E[ei(pi)7y?] = a2(2E[ei(pi)77i] + E[e'i(pi)r;?]). (69) 

Thus, ( l66b reads 

E[ei(pi)^i?7^] =a,ala^E[e[{p,)7i] + E[ei(pi)77i] 

X (aiE[^] - aiaja^ + 2a2E[?ji^2]) (70) 
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(71) 



which, by using (|9]l, can also be reexpressed as 

In turn, we have 

E[eUpi)^] =a^E[eUpi)?7?] +2a2E[el(pi)?7i]Efe] 
+ E[e;(pi)]E[r)i] 

=«^E[eUpi)^?] + E[e;(pi)](E[^] ~ ala^) 

(72) 

which, by using ( |57] l, leads to 

E[e'i(pi)?^]E[r;i?;i] ^aia^^^^jQ, (^^^^2] 

+ a2E[el(pi)]ai(E[^yf]-a2^2N 



(73) 



From the difference of dTTT l and ( l73T l. we derive that 



which, by using again (|57] l, yields (fT2l l. 

Finally, we will prove Formula (fTsT l. We decompose 771 as 
follows: 

Vi ^ br]2 +V1 where = £[^1^2], (75) 

cr^ — E[?7|] and, 77!^ is independent of (?72, ui, ^2)- This 
allows us to write 

mi{pi)Q2iP2)mm] = b E[ei{pi)Q2{p2)vl] 

+ E[eMe2iP2)vtm]- (76) 

Let us first calculate E[9i(pi)9(/92)'72]- For * ^ {1)2}, 
consider the decomposition: 



(77) 



where CiCr^ = E[?7i?72] and, ry2, (Vi^vi^) ™d {vi,V2) are 
independent. We have then 

E[ei{pi)e2{p2)vi I '7j^,»7^,''^i,f^2] 
1 



ei(ui + ciC + r/i )02(t^2 + caC + vt) 



C'exp(-^)dC. (78) 



and, for every (ti, r2) € R^, lim|^|^oo 0i(n + ciC)02(t-2 + 
C2C)C6xp ( — ^5-) = since, for |C| large enough. 



"Icll&liTi + CiC)Q2{t2 + C2C)CI exp ( - 



< |ei(ri+ciC)|(ciC)'exp(-^) 

X |e2(T2 +C2C)|(c2C)^exp I 



(80) 



and Assumption [(i)] holds. We can therefore deduce, by in- 
tegrating by parts in (iTSl l and taking the expectation w.r.t. 

(77i",'?2",f^i,W2), that 

E[ei(pi)e2(p2)»?^] 

=3?2(E[ei(pi)e2(p2)] + ciE[e'i(pi)e2(p2)^2] 

+ c2E[ei(pi)e^(P2)?72]) 
=52E[ei(pi)e2(p2)] + E[e'i(pi)e2(p2)72]E[77i??2] 

+ E[ei(pi)e'2(p2)72]E[r72772]. (81) 

Let us now calculate E[8i(pi)82(p2)'7i"'72]- We have, for i e 
{1, 2}, rjt = Ctrji + fjl, where 



(82) 



It can be noticed that 



E[|ei(pi)e2(p2) + (cie'i(pi)e2(p2) + c2Qi{pi)Q'2{p2))m\] 

< E.[\QM\Y'^E.[\Q2{P2)?Y" + |ci|(E[|ei(pi)|^]^/^ 

X E[|e2(p2)h'']^/''' + |c2|E[|ei(pi)p]V3E[|e^(P2))p]^/^) 

xE[|72p]'/'<oo (79) 



and, 77^'- is independent of (772, iji ,7^2 ,vi,V2)- By proceeding 
similarly to the proof of dsTl ). we get 

E[ei(pi)e2(p2)^7^ I m,fii,fi^] = E[i?iif] 
X (ciE[e;(pi)e2(p2) 1^72,^1^,%^] 

+ c2E[ei(pi)ej,(p2) I m,fit,V2]) (83) 

which, owing to ( |82] |. allows us to write 

E[ei(pi)e2(p2)?y^72] =E[e;(pi)e2(P2)572]E[r7i?7j^] 

+ E[ei(pi)ei(p2)??2]E[ry277^]. 

(84) 

On the other hand, from ( fTSl ), we deduce that, for i e {1,2}, 
^[ViVi] = E[7/i7i] - 6E[?7j72], so yielding 

E[ei(pi)e2(p2)^?^72] = E[e'i(pi)e2(p2)72]E[77i,7i] 

+ E[ei(pi)ei(p2)72]E[r72m] 
^6(E[e;(pi)e2(p2)72]E[r;i572] 

+ E[ei(pi)e^(p2)72]E[r7272]). 

(85) 

Altogether, ^ and ^ lead to 

E[6i(pi)92(p2)??i72] = E[6i(pi)62(p2)]E[7if72] 

+ E[e;(pi)e2(p2)72]E[r/i7i] + E[ei(pi)eJ,(p2)?72]E[r/2?7i]- 

(86) 

In order to obtain a more symmetric expression, let us 
now look at the difference E[8i(pi)89(p2)?72]E[7727i] ~ 
E[ei(pi)e^(p^2)7i]E[?725?22 = E[ei(pi)e^,(p2)57i2] where 
7712 — E[77277i]?72 — E[772?72]77i. Since 7712 is a linear combination 
of rji and 772 and, E [772 771 2] = 0, 7712 is independent of 
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{'i]2,vi,V2)- Similarly to the derivation of Formula (fTOl i. it 
can be deduced that 

E[ei(pi)e2(p2)?7l2 I V2,Vl,V2] 

=£[61(^1)7712 1 ?72, ui, W2] 82(^2 + m) 

=E[ei(pi) I ?72,ui,t^2]E[77i 7712] 62(^2) (87) 
which leads to 

E[ei(pi)e^(p2)572]E[r;2^i] - E[ei(pi)e^(p2)^7i]E[r;2?72] 

=E[ei(pi)ei(P2)^7i2] - E[e;(pi)e'2(p2)]E[??i^7i2] 

=E[e'i(pi)e^(p2)](E[??i?72]E[r?2^?i] - E[77i57i]E[7y2^2]). (88) 
Eq. (fTsT l is then derived by combining ( |86] | with dSST l. 

Appendix B 
Proof of Proposition!!] 
We have, for every p G D, 

E[|i?(p)p] = E[|5(p)p] + E[|iV(p)p]+2Re{E[iV(p)(5(p))*]} 

(89) 

Since n and s are uncorrelated, this yields 
E[|5(p)|2] = E[|i?(p)|2] - E[|A^(p)|2]. In^ addition, we have 
E[5(p)(5(p))*] ^ E[i?(p)(5(p))*] - E[7V(p)(5(p))*]. The 
previous two equations show that 

Vp e B, E[|^(p) - S{p)f] = Ejl^(p) - i?(p)n 

-E[|7V(p) p] + 2 Re{E[Ar(p) (§(p)) *]}. (90) 

Moreover, using ( fTTI l, the second term in the right-hand side 
of dinil is 



Combining this equation with ( |92] ) and ( |94l ) yields 

Emp)iSip)y]^,±Ei&,ir,) f^^^^^^ (96) 



«=i 



Gathering now ( |90l l, ( |9T] ) and ( |96] l, ( fTSl l is obtained. 

From Parseval's formula, the global MSE can be expressed 

as 

D E[£{s- s)] = ^Yl E[I^(P) - S{p)n (97) 



D 
peD 

The above equation together with ( fTST i show that ( fT9] l holds 
with 



DA = 7(2^ E[e^(r,)] Rc{7j - l^(P)r') • (98) 

Furthermore, by defining 

V£e {!,...,£}, = (99) 

and using ( fT6b . it can be noticed that 

E[nfn£]= ^ E[n(y)n(x)]^f (y)¥3f(x) 

(x,y)eD2 

^ E[iV(p)(7V(p'))*] 



1)2 



(P,P')6I 

pel 



7 ^ $^(p)($£(p))' 



D ^ 



77£- 



(100) 



E[|iV(p)|^ 



(91) 



|i7(p)p- 

On the other hand, according to (O, the last term in the right- 
hand side of (l90l l is such that 

L 

E[7V(p)(5(p))*]=EEE[^^«(x)] 

xGD £=1 

X exp(-27r7x^D~V)(5«(p))*- (92) 

Furthermore, we know from ^ that = + ?^^), where 

ni = (n, (^^), ui, = {u, (fe) (93) 

and, u is the field in R-Di x - xr>d whose discrete Fourier 
coefficients are given by ((Sj. From (|93] l as well as the 
assumptions made on the noise n corrupting the data, it is 
clear that (n^, n(x)) is a zero-mean Gaussian vector which is 
independent of ui. Thus, by using (fTOl i in Proposition |2] we 
obtain: 

E[?,H(x)] = E[eKr,)]E[n,H(x)]. (94) 
Let us now calculate E[nf n(x)]. Using (|93] | and ( fTSI l. we get 

E[n£n(x)] = E E[n(x)r7,(y)]93^(y) 

yGD 



Hence, as claimed in the last part of our statements, {ji)i<i<L 
is real-valued since it is the cross-correlation of a real-valued 
sequence. 

Appendix C 
Proof of Proposition^ 

By using (|25]l, we have X^xgd^W^W = J^i^i^ene, 
where has been here redefined as 



ni = {n,ip). 



(101) 



In addition, £{n) = i?"" Epeg l^(p)l Vl^(p)l'- This al- 
lows us to rewrite ( [34l i as £0 — £0 ^ —2A — B + 2C, where 



(102) 



pen; 



J_ \N{p)l/D^ 



P6 

L 



(103) 



(104) 



^ E 
(p'.p'Oen: 



E [iV(p') (7V(p")) *] exp(2^ix^r>-ip') 



^pt;^(p') 



(95) 



The variance of the error in the estimation of the risk is thus 
given by 

yar[£o-£o] = E[(£o-£o)2] =4E[A2]+4E[AB]-8E[AC7] 
+ E[b2]-4E[BC]+4E[C^]. (105) 
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We will now calculate each of the terms in the right hand-side 
term of the above expression to determine the variance. 

• Due to the independence of s and n, the first term to 
calculate is equal to 



[^'] = 7^ E E[5(p)(5(p'))lE[iV(p)(7V(p'))*: 

(106) 



(p,p')eQ" 

By using ST% . this expression simplifies as 



7 ^ E[|5(p)p] 



E 



i?3^ |i7(p)|2 



(107) 



The second term cancels. Indeed, since n and hence n 
are zero-mean, 

^[^^1 = 4 E E[.(x)]E[n(x)(n(x'))^]. (108) 



(x,x')ei 



and, since (n(x),n(x')) is zero-mean Gaussian, it has a 
symmetric distribution and E[n(x)(n(x'))^] — 0. 
The calculation of the third term is a bit more involved. 
We have 

1 ^ 

E[^C'] = -j^^(E[s(x)5,n,n(x)] 



-E[,s(x)eKr,)n(x)]77,). (109) 



In order to find a tractable expression of 
E[s(x)s^ n(x)] with £ E {l,...,L}, we will 
first consider the following conditional expectation w.rt. 
s: E[s{x)'senin{-x.) \ s] — s(x) E[8(r^)ri£ H(x) | s]. 
According to Formula dTTT) in Proposition l2n ~ 

E[0{ri)nen{x.) \ s] =E[e^(rf )n(x) | s] E[neni] 

+ E[ee{n) I s]E[nen{K)] (110) 

which, by using dlOOt , allows us to deduce that 

E[s(x)e(r,)n,n(x)] = E[s(x)eKr,)n(x)] 77, 

+ E[s(x)e,(r,)]E[H,n(x)]. (Ill) 



This shows that ( |109t can be simplified as follows: 

1 ^ 

E[^^] = ;^EEE[Kx)e,(r,)]E[n,n(x)]. (112) 
Furthermore, according to ( fTTj l and ( llOlt . we have 



E[n,n(x)] =-i^ J2 E[7V(p)(iV(p'))1($.(p))' 

(P,P')GQ2 

X exp(-27rix^D^V) 

z 

D 



7 (^^(P))* . „ Tn-i N 
iu/„M2 exp(-27rix D p) 



pet; 



l^(p) 



(113) 



*Proposition|2]is applicable to the calculation of the conditional expectation 
since conditioning w.r.t. s amounts to fixing ue (see the remark at the end of 
Section inn. 



This yields 

E\AC] 



peQ e=i 



—y 

peQ 



_ \Hip) 
E[5(p)(5(p))*] 



I^(P) 



(114) 



The calculation of the fourth term is more classical since 
\N{p)\'^ / D is the p bin of the periodogram 115111 of the 
Gaussian white noise n. More precisely, since \N{p)\'^ / D 
is an unbiased estimate of 7, 



E[B^ 



1 



E 

(p,p')e( 



Cov(|iV(p)|M7V(p')| 



|i/(p)P|i/(p 



(115) 



In the above summation, we know that, if p 7^ p' and p 7^ 
Dl - p' with 1 = (1, . . . , 1)^ e R'^, N{p) and N{p') 
are independent and thus, Cov(|7V(p)p, |iV(p')p) = 0. 
On the other hand, if p = p' or p — Dl — p', then 

Cov(|iV(p)|2, |iV(p')P) = E[|iV(p)|4] - ^^D\ Let 

S = {p = (pi,---,Pd)^ee| 

Vze {!,..., d},p,e {0, A/2}}. (116) 

If p e §, then A^(p) is a zero-mean Gaussian real random 
variable and E[|A^(p)|4] = 3E[{N{p))Y = ^l^D^. 
Otherwise, A^(p) a zero-mean Gaussian circular com- 
plex random variable and E[|iV(p)|''] = 2E[| A^(p)|2]2 
= It can be deduced that 



E[i?2] 



E 

, peQns 



E ( 

pGQn(D\S) 



Var[(iV(p))' 



I^(P)I^ 

Var[|iV(p)|^] 
I^(P)I^ 
Cov(|jV(p)|MjV(m-p)p) 

|F(p)P|i7(m-p)p 
^ / 272152 



> pGQnS 



I^(P) 



E ( 



pGQn(D\S) 



72^2 
H{P)\ 



72/^2 

\h{pW 



2f_ 

£12 



E 



1 



pen; 



^(P)l 



4 ■ 



(117) 



Let us now turn our attention to the fifth term. According 
to ( [Tol l and the definition of 7^ in ( llOOI l, for every I e 
{1, . . . , L}, s'e rie — Q'i{re)'yji is zero-mean and we have 
then 



xGD e=i 



E[e^(r,)(n(x))']77,). (118) 
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By applying now Formula ST% in Proposition |2l we have, 
for every £ G {1, . . . , L}, 

E[5,n,(n(x)) V E[eKr,)(n(x))']77, 
=E[e,(r£)n,(H(x))'] - E[e^(r,)(n(x))'] E[n,n,] 
=2E[e',{re)] E[n,n(x)] E[n(x)n,] (119) 

where, in compliance with (|24] |. Ui is now given by 

ne^(n,ipi). (120) 

Furthermore, similarly to ( |95l ), we have 

E[n(x)n,]=ij E E[7V(p)(7V(p'))1$.(p') 



(p,p')et 



X exp(27rix £) p) 

P'6Q ' 



Altogether, (fTTlT l, (fTT9] l and (fTITl l yield 

E[S-,n,(n(x))'] - E[eKr,)(n(x))']77, 

27V.^,, ^ 'i>^(pO('i>^(p))^ 

= ;^E[e,(r,)] ^ 



(121) 



i/(p')|i?(p)P 
X exp {2mx^D-\p' - p)). (122) 



Hence, jl 18l l can be reexpressed as 

27^ ' 



£=1 



where 



_ 1 ^ $,(p)(<i>,(p))* 
''^ i?(p)|i^(p)P 

Let us now consider the last term 

L L 

D2" 



(124) 



^ L L 

^t*^^] 752 X! X! (E[s£Sin£n,] - E[Sie^(r£) n^]77f 



« = 1 2=1 



E[s£eKrOn,]772 + E[e^(r,)eKrO]7'm) • (125) 



Appealing to Formula ( fTSl l in Proposition |2] and (llOOI l. 
we have 

E[e^(r^)ei(ri)n£ni] 

+ E[ee{n)e',{n)nehi, + E[e',{n)e'M)] 

X (E[nine\E[rieni\ - 7^7<>7j). (126) 

This allows us to simplify ( 11251 ) as follows: 



L L 



E[^']=;^EE(E[^^?.]E[n£«.] 

1=1 1=1 

+ E[Q',{r,)Q[{n)]E[n,ni\E[riin,]). (127) 



Furthermore, according to ( fTTj i, dlQll i, ( fT6] l and ( 11201 ), we 
have 

= E E[iV(p)(iV(p'))1(5.(p))*$.(p') 
(p,p')eQ' 

($,(p))*$,(p) 



E 



PG4 



l^(p)l^ 



(128) 



and 



EK^d-;^ E e[7V(p)(7V(p'))>Kp')(<I'.(p))* 



(p,p')et 



(129) 



where the expression of 7^ ^ is given by dJTl ). Hence, by 
using ( ll28l )-( fT29] l, ( I127l i can be rewritten as 



pel 

„2 



Ii/(P)P 

L L 



EE[0£(^^)0U^O]7i.7.,^- (130) 



= 1 i=l 



• In conclusion, we deduce from ( 11051 ). (1107b . (1114b . (1117b . 
( fT23T l and ( fTIOl l that 

Var[^o-.^.]^i;iE''''^"^~'^^^'^^ 



§(E EE[©K^^)0Kn)]7, 



1)3 

pen 



|i?(p) 



-2EE[eM.. + ^E^)- (131) 

£=1 peQ ' ^^'^ 

By exploiting now (fTsT l (see Proposition |4|i and noticing 
that (k£)i<£<l is real-valued, this expression can be 
simplified as follows: 



|i7(p)P 



472 



pet 

L L 



-1 1=1 



27 V 1 

^ li/(p)|4- 



(132) 



pet; 



Eq. ( l35T l follows by using Parseval's formula. 
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